
tstart=1;
if length(incomevals)==3
tfinish=25;
elseif length(incomevals)==2
tfinish=10;
end
weight=0;
pen=1;
if beta==0.86
    pen=4;
elseif beta==0.88
    pen=3;
elseif beta==0.9 
    pen=3;
elseif  beta==0.92
    pen=2;
end  
%beta=0.86 needs a penalty of 0.04 to achieve stability, beta=0.88 needs a penalty of 0.03 to achieve stability%
%beta=0.90 needs a penalty of 0.03 to achieve stability, beta=0.92 needs a penalty of 0.02 to achieve stability%
%beta=0.94 and beta=0.96 does not need a penalty;


if Qbeta==3 & length(incomevals)==2
pen=1; 
elseif Qbeta>=3 & length(incomevals)==2 %%%this is the penalty for beta=0.9 when there are only 2 income states
pen=1; %%%this is the penalty for beta=0.92 and beta=0.94 and beta=0.96 when there are only 2 income states 
end
global gammagrid elementsize lowerbound table m n agent
agent=3;
beta=betavec(Qbeta);
if beta<=0.93
elementsize=0.1;
lowerbound=0.1;
[gammagrid,m]=grid3(elementsize,lowerbound);
else
elementsize=0.05;
lowerbound=0.05;
[gammagrid,m]=grid3(elementsize,lowerbound);
end

gammagrid(:,4)=1;
gammagrid(:,5)=gammagrid(:,2)./gammagrid(:,1);
gammagrid(:,6)=gammagrid(:,3)./gammagrid(:,1);
n=length(gammagrid);
table=[]
for r1=1:m
index1=sum(m:-1:m-(r1-2));
for r2=1:m-(r1-1)
table(r1,r2)=index1+r2;
end
end

d1=n;   %length of gammagrid
d2=length(incomevals)^3;  %states
d3=3;   %agent


WSB=ones(d1,d2,d3)*1000;
PHI=ones(d1,d2,d3)*1000;
CONS=ones(d1,d2,d3)*1000;

XXKEEP=ones(7,d2,d3)*1000;
GAMMAGRID=gammagrid;
PROB=ones(d2,d2)*1000;
STATE=ones(d2,d3)*1000;
TABLE=table;
COUNTGRID=ones(d1,d2);
ADMIT=ones(7,d2); %%%this is the matrix where you check whether it is an admissible deviation



income=incomevals;

global P beta S state agent Y scalefactor rho ewdev statedev12 statedev13 statedev23 wdev cdev
beta=betavec(Qbeta);
PEN=[0.1 0.15 0.5 1];
penalty=0+0.01*(pen-1);
S=gridmake(income, income, income);
SS=gridmake([1:length(income(:,1)) ]', [1:length(income(:,1)) ]', [1:length(income(:,1)) ]');

state=length(S);
Y=sum(S,2);				% possible aggregate income outcome

for k=1:state
    for j=1:state
        P(j,k)=Q1(SS(j,1),SS(k,1))*Q1(SS(j,2),SS(k,2))*Q1(SS(j,3),SS(k,3));
    end
end
agent=3;  
caut=zeros(state,agent);
% Compute set of autarky values
for k=1:agent
caut(:,k)=S(:,k);
end
waut=zeros(state,agent);
wautnew=waut;
dwaut=1;
crit=0.0000000001;
while dwaut>crit
    
    wautnew(:,1)=utility(rho,caut(:,1))+beta*P*waut(:,1);
    wautnew(:,2)=utility(rho,caut(:,2))+beta*P*waut(:,2);
    wautnew(:,3)=utility(rho,caut(:,3))+beta*P*waut(:,3);
    dwaut=max(max(abs(wautnew-waut)));
    waut=wautnew
end
for j=1:state
    for g=1:agent
    if waut(j,g)>0
waut(j,g)=waut(j,g)*(1-penalty);
    else
waut(j,g)=waut(j,g)/(1-penalty);
    end
    end
end

w=zeros(n,state,agent);
gammar=zeros(n,state,agent);
c=zeros(n,state,agent);

% Note that phi contains the Lagrange multipliers on the enforcement
% constraint
phi=zeros(n,state,agent);

count=0

for j=1:state
    for g=1:agent
c(:,j,g)=Y(j)*gammagrid(:,agent+g)./(1+gammagrid(:,5)+gammagrid(:,6));
     end
end
cfirstbest=c;

gammar=zeros(n,state,agent);
for j=1:state
for k=1:agent   
gammar(:,j,k)=gammagrid(:,agent+k);

end
end
% Compute the value functions
w=zeros(n,state,agent);
wnew=zeros(n,state,agent);
dw=1;
while dw>crit
    for j=1:state
        wnew(:,j,1)=utility(rho,cfirstbest(:,j,1))+beta*w(:,:,1)*P(j,:)';
        wnew(:,j,2)=utility(rho,cfirstbest(:,j,2))+beta*w(:,:,2)*P(j,:)';
        wnew(:,j,3)=utility(rho,cfirstbest(:,j,3))+beta*w(:,:,3)*P(j,:)';

    end
            dw=max(max(max(abs(w-wnew))));
        w=wnew;
end
wfirstbest=w;
%%%
ewfirstbest=zeros(n,state,agent);
for j=1:state
        for g=1:agent
ewfirstbest(:,j,g)=wfirstbest(:,:,g)*P(j,:)';
        end
end
ewsb=ewfirstbest;
wsb=wfirstbest;
wsbnew=wsb;

% Now we can start iterating to update our guesses of the value functions
% when there are enforcement constraints. 
%%%%Load in the results for group of 2
load('RESULTS2ALL','WDEV2','PROBDEV2','STATEDEV2','EWDEV2');
wdev=real(WDEV2(:,:,:,Qbeta,Qrho));
Pdev=PROBDEV2(:,:,1);
Sdev=STATEDEV2(:,:,1);
ewdev=real(EWDEV2(:,:,:,Qbeta,Qrho));
cdev=real(CDEV2(:,:,:,Qbeta));
%clear WDEV2 PROBDEV2 STATEDEV2 EWDEV2
statedev=size(Sdev,1);
idev=size(wdev,1);
agentdev=2;
for i=1:idev
for j=1:statedev
    for g=1:agentdev
    if ewdev(i,j,g)>0
ewdev(i,j,g)=ewdev(i,j,g)*(1-penalty);
wdev(i,j,g)=wdev(i,j,g)*(1-penalty);
    else
ewdev(i,j,g)=ewdev(i,j,g)/(1-penalty);
wdev(i,j,g)=wdev(i,j,g)/(1-penalty);
    end
    end
end
end
%find the respective states for 2 people
for k=1:state
    statedev12(k,1)=find((Sdev(:,1))<=(S(k,1))+0.00001 & (Sdev(:,1))>=(S(k,1))-0.00001 & (Sdev(:,2))>=(S(k,2))-0.00001 & Sdev(:,2)<=S(k,2)+0.00001);
    statedev13(k,1)=find((Sdev(:,1))<=(S(k,1))+0.00001 & (Sdev(:,1))>=(S(k,1))-0.00001 & (Sdev(:,2))>=(S(k,3))-0.00001 & Sdev(:,2)<=S(k,3)+0.00001);
    statedev23(k,1)=find((Sdev(:,1))<=(S(k,2))+0.00001 & (Sdev(:,1))>=(S(k,2))-0.00001 & (Sdev(:,2))>=(S(k,3))-0.00001 & Sdev(:,2)<=S(k,3)+0.00001);
end


gap=1;
gapend=0.001;
options=optimset('Display','iter');
scalefactor=1;
%scalefactor=20/range(ewfirstbest(:,2,1));
count=0;
eps=0.0000001;
t=0;
load('RESULTS3AUT','CONSAUT','XXKEEPAUT');

agentperm=perms(1:agent);
history=zeros(n,state);
gridperm=[]; iperm=[]; jperm=[]; csol=[]; wsbnewsol=[]; phisol=[]; outputkeep=[]; xkeep=[]; xkeep13=[]; outputkeep13=[]; xkeep23=[]; outputkeep23=[]; xkeepperm=[]; outputkeepperm=[];
for t=1: tfinish
countgridod=ones(length(ewdev),27);
for j=1:state
if one_per_aut==1
Optimaldeviation_perm_2018_05_10;
else
Optimaldeviation_perm_2019_09_09;
end
end
agentperm=perms(1:agent);
count=0
bla1=zeros(n,state);
bla2=zeros(n,state);
countgrid=[ones(n,state)];

for j=1:state
stateperm=perms(S(j,:));

%find the admissible deviations
%%%%consider deviations of 1 and 2, we want to check if 13 or 23 have a
%%%%better deviation. Keep 1 same, see if 3 wants to swap in, Keep 2 same,
%%%%see if 3 wants to swap in. 
admissible=[];
for d=1:idev
    A=ones(idev,1)*[outputkeep(d,j,1)  outputkeep(d,j,3)];
    B=[outputkeep13(1:idev,j,1) outputkeep13(1:idev,j,3)];
    C=sum(A+eps<=B,2);
    rr13=find(C<2);
    A=ones(idev,1)*[outputkeep(d,j,2)  outputkeep(d,j,3)];
    B=[outputkeep23(1:idev,j,2) outputkeep23(1:idev,j,3)];
    C=sum(A+eps<=B,2);
    rr23=find(C<2);
   if size(rr13,1)==idev & size(rr23,1)==idev
       admissible=[admissible; d];
   end
end
    
admissible13=[];
for d=1:idev
    A=ones(idev,1)*[outputkeep13(d,j,1)  outputkeep13(d,j,2)];
    B=[outputkeep(:,j,1) outputkeep(:,j,2)];
    C=sum(A+eps<=B,2);
    rr12=find(C<2);
    A=ones(idev,1)*[outputkeep13(d,j,2)  outputkeep13(d,j,3)];
    B=[outputkeep23(:,j,2) outputkeep23(:,j,3)];
    C=sum(A+eps<=B,2);
    rr23=find(C<2);
   if size(rr12,1)==idev & size(rr23,1)==idev
       admissible13=[admissible13; d];
   end
end

admissible23=[];
for d=1:idev
    A=ones(idev,1)*[outputkeep23(d,j,1)  outputkeep23(d,j,2)];
    B=[outputkeep(:,j,1) outputkeep(:,j,2)];
    C=sum(A+eps<=B,2);
    rr12=find(C<2);
    A=ones(idev,1)*[outputkeep23(d,j,1)  outputkeep13(d,j,3)];
    B=[outputkeep13(:,j,1) outputkeep13(:,j,3)];
    C=sum(A+eps<=B,2);
    rr13=find(C<2);
   if size(rr12,1)==idev & size(rr13,1)==idev
       admissible23=[admissible23; d];
   end
end
admissible, admissible13, admissible23
for i=1:n
if countgrid(i,j)>0    
index=[i,j];

for d=1:idev
   gamma2implied_12(d,1)=CDEV2(d,statedev12(j),1,Qbeta)/(S(j,1)+S(j,2));
   gamma2implied_13(d,1)=CDEV2(d,statedev13(j),1,Qbeta)/(S(j,1)+S(j,3)); 
   gamma2implied_23(d,1)=CDEV2(d,statedev23(j),1,Qbeta)/(S(j,2)+S(j,3)); 
end
gamma2implied_12(:,2)=1-gamma2implied_12(:,1);
gamma2implied_13(:,2)=1-gamma2implied_13(:,1);
gamma2implied_23(:,2)=1-gamma2implied_23(:,1);
%%%%find joint deviations that are possible

dev12=sum(ones(idev,1)*[squeeze(utility(rho,cfirstbest(i,j,1:2))+beta*ewsb(i,j,1:2))'] -squeeze(wdev(:,statedev12(j),1:2))<0,2);
dev13=sum(ones(idev,1)*[squeeze(utility(rho,cfirstbest(i,j,1:2:3))+beta*ewsb(i,j,1:2:3))']-squeeze(wdev(:,statedev13(j),1:2))<0,2);
dev23=sum(ones(idev,1)*[squeeze(utility(rho,cfirstbest(i,j,2:3))+beta*ewsb(i,j,2:3))']-squeeze(wdev(:,statedev23(j),1:2))<0,2);

coalitionbind12=find(dev12==2);
coalitionbind13=find(dev13==2);
coalitionbind23=find(dev23==2);

gridperm=perms(gammagrid(i,1:3));
for k=1:length(gridperm)
iperm(k)=find(gammagrid(:,1)<=gridperm(k,1)+eps & gammagrid(:,1)>=gridperm(k,1)-eps& gammagrid(:,2)<=gridperm(k,2)+eps & gammagrid(:,2)>=gridperm(k,2)-eps & gammagrid(:,3)<=gridperm(k,3)+eps & gammagrid(:,3)>=gridperm(k,3)-eps);
jperm(k)=find(S(:,1)==stateperm(k,1) & S(:,2)==stateperm(k,2) & S(:,3)==stateperm(k,3));
end
% Constraint binding for
% Person 1 only

if scalefactor*(utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))<scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
        & scalefactor*(utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))>=scalefactor*(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
        & scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
        & dev12<2  & dev13<2 & dev23<2
   
%%%%check whether it's possible that the individual deviation goes through    
ccheck(2)=(Y(j)-xkeep(1,j,1))/(1+gammagrid(i,3)/gammagrid(i,2));
ccheck(3)=Y(j)-xkeep(1,j,1)-ccheck(2);

if ccheck(2)>=0.9*xkeep(1,j,2) & ccheck(3)>=0.9*xkeep13(1,j,3)
%%%auxiliary problem to get good starting guess
clear x0
xaux2 = (Y(j)-xkeep(1,j,1))/(1+gammagrid(i,6)/gammagrid(i,5));
xaux3 = gammagrid(i,6)/gammagrid(i,5)*xaux2;
x0(1)=Y(j)-xaux2-xaux3;
x0(2)=xaux2;

x0=[x0 x0(2)/x0(1) (Y(j)-x0(2)-x0(1))/x0(1)];
        
[xsol, output]=agent1solverho(x0,index,ewsb,waut);
if output(6)<=0
   error('did not converge') 
end
csol(i,j,1:2)=xsol(1:2);
csol(i,j,3)=Y(j)-xsol(1)-xsol(2);
gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
wsbnewsol(i,j,1:agent)=output(1:agent);
phisol(i,j,1)=(gammagrid(i,5)/gammar(i,j,2))^rho-1;
phisol(i,j,2:agent)=0;
%%%%%here you then need to check that it goes through ex-post!
if csol(2)<xkeep(1,j,2) 
if length(admissible)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep(admissible,j,:)));
else
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep(admissible,j,:))');
end

history(i,j)=admissible(d);
csol(i,j,1)=xkeep(admissible(d),j,1);
csol(i,j,2)=xkeep(admissible(d),j,2);
wsbnewsol(i,j,:)=outputkeep(admissible(d),j,1:agent);
if equalsharing==1
history(i,j)=4;
csol(i,j,1)=xkeep(4,j,1);
csol(i,j,2)=xkeep(4,j,2);
wsbnewsol(i,j,:)=outputkeep(4,j,1:agent);
end   
phisol(i,j,3)=0;


csol(i,j,3)=Y(j)-csol(i,j,1)-csol(i,j,2);
gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
phisol(i,j,1)=(gammagrid(i,6)/gammar(i,j,3))^rho-1;
phisol(i,j,2)=((gammagrid(i,6)/gammagrid(i,5))/(gammar(i,j,3)/gammar(i,j,2)))^rho-1;

elseif csol(3)<xkeep13(1,j,3)
if size(admissible13,1)>1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep13(admissible13,j,:))');
else
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep13(admissible13,j,:)));
end
csol(i,j,1:2)=xkeep13(admissible13(d),j,1:2);
csol(i,j,3)=Y(j)-xkeep13(admissible13(d),j,1)-xkeep13(admissible13(d),j,2);
wsbnewsol(i,j,1:agent)=outputkeep13(admissible13(d),j,1:agent);
if equalsharing==1
csol(i,j,1:2)=xkeep13(4,j,1:2);
csol(i,j,3)=Y(j)-xkeep13(4,j,1)-xkeep13(4,j,2);   
wsbnewsol(i,j,1:agent)=outputkeep13(4,j,1:agent);
end
gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);

phisol(i,j,1)=(gammagrid(i,5)/gammar(i,j,2))^rho-1;
phisol(i,j,3)=((gammar(i,j,3)/gammar(i,j,2))/(gammagrid(i,6)/gammagrid(i,5)))^rho-1;
phisol(i,j,2)=0;    
end
%%%%here you check that it goes through ex-ante 
elseif ccheck(2)<0.9*xkeep(1,j,2) 
if length(admissible)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep(admissible,j,:)));
else
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep(admissible,j,:))');
end
history(i,j)=admissible(d);
csol(i,j,1)=xkeep(admissible(d),j,1);
csol(i,j,2)=xkeep(admissible(d),j,2);
wsbnewsol(i,j,:)=outputkeep(admissible(d),j,1:agent);
if equalsharing==1
history(i,j)=4;
csol(i,j,1)=xkeep(4,j,1);
csol(i,j,2)=xkeep(4,j,2);
wsbnewsol(i,j,:)=outputkeep(4,j,1:agent);   
    
end
phisol(i,j,3)=0;

csol(i,j,3)=Y(j)-csol(i,j,1)-csol(i,j,2);
gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
phisol(i,j,1)=(gammagrid(i,6)/gammar(i,j,3))^rho-1;
phisol(i,j,2)=((gammagrid(i,6)/gammagrid(i,5))/(gammar(i,j,3)/gammar(i,j,2)))^rho-1;

elseif ccheck(3)<0.9*xkeep13(1,j,3)
if size(admissible13,1)>1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep13(admissible13,j,:))');
else
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep13(admissible13,j,:)));
end
csol(i,j,1:2)=xkeep13(admissible13(d),j,1:2);
csol(i,j,3)=Y(j)-xkeep13(admissible13(d),j,1)-xkeep13(admissible13(d),j,2);
wsbnewsol(i,j,1:agent)=outputkeep13(admissible13(d),j,1:agent);
if equalsharing==1
csol(i,j,1:2)=xkeep13(4,j,1:2);
csol(i,j,3)=Y(j)-xkeep13(4,j,1)-xkeep13(4,j,2);   
wsbnewsol(i,j,1:agent)=outputkeep13(4,j,1:agent);
end

gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);

phisol(i,j,1)=(gammagrid(i,5)/gammar(i,j,2))^rho-1;
phisol(i,j,3)=((gammar(i,j,3)/gammar(i,j,2))/(gammagrid(i,6)/gammagrid(i,5)))^rho-1;
phisol(i,j,2)=0;
elseif ccheck(2)<1*xkeep(1,j,2)  & ccheck(3)<1*xkeep13(1,j,3) 
 
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep(admissible,j,:))');
history(i,j)=admissible(d);
csol(i,j,1)=xkeep(admissible(d),j,1);
csol(i,j,2)=xkeep(admissible(d),j,2);
wsbnewsol(i,j,:)=outputkeep(admissible(d),j,1:agent);
if equalsharing==1
history(i,j)=4;
csol(i,j,1)=xkeep(4,j,1);
csol(i,j,2)=xkeep(4,j,2);
wsbnewsol(i,j,:)=outputkeep(4,j,1:agent);   
  
end
phisol(i,j,3)=0;
bla1(i,j)=1;
csol(i,j,3)=Y(j)-csol(i,j,1)-csol(i,j,2);
gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
phisol(i,j,1)=(gammagrid(i,6)/gammar(i,j,3))^rho-1;
phisol(i,j,2)=((gammagrid(i,6)/gammagrid(i,5))/(gammar(i,j,3)/gammar(i,j,2)))^rho-1;

end

countgrid(i,j)=0;
% identical cases
for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
countgrid(iperm(k),jperm(k))=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constraint binding for Person 1 and Person 2%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif  size(coalitionbind12)>0;

if length(admissible)==1
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep(admissible,j,:)));
else   
[CC,d]=max([(gammagrid(i,1)-0.0001) gammagrid(i,2:3)]*squeeze(outputkeep(admissible,j,:))');
end
history(i,j)=admissible(d);
csol(i,j,1)=xkeep(admissible(d),j,1);
csol(i,j,2)=xkeep(admissible(d),j,2);
wsbnewsol(i,j,:)=outputkeep(admissible(d),j,1:agent);
if equalsharing==1
history(i,j)=4;
csol(i,j,1)=xkeep(4,j,1);
csol(i,j,2)=xkeep(4,j,2);
wsbnewsol(i,j,:)=outputkeep(4,j,1:agent);   
  
end
phisol(i,j,3)=0;

csol(i,j,3)=Y(j)-csol(i,j,1)-csol(i,j,2);
gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
phisol(i,j,1)=(gammagrid(i,6)/gammar(i,j,3))^rho-1;
phisol(i,j,2)=((gammagrid(i,6)/gammagrid(i,5))/(gammar(i,j,3)/gammar(i,j,2)))^rho-1;

countgrid(i,j)=0;
% identical cases
for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
countgrid(iperm(k),jperm(k))=0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Constraint binding for noone%%%%%%%
elseif scalefactor*(utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))>=scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
        & scalefactor*(utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))>=scalefactor*(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
        & scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
        & dev12<2  & dev13<2 & dev23<2
csol(i,j,:)=cfirstbest(i,j,:);
gammar(i,j,:)=gammagrid(i,4:6);
phisol(i,j,1:agent)=0;
for g=1:agent
wsbnewsol(i,j,g)=utility(rho,cfirstbest(i,j,g))+beta*ewsb(i,j,g);
end
countgrid(i,j)=0;

% identical cases
for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
countgrid(iperm(k),jperm(k))=0;
end


end
end % end of countgrid if loop
end % end of i loop
end % end of j loop



for g=1:agent
gap(g)=scalefactor*norm((wsbnew(:,:,g)-wsb(:,:,g)));
end
for j=1:state
for g=1:agent
gammar(:,j,g)=c(:,j,g)./c(:,j,1);
end
end
gap=sum(gap)
wsb=wsbnew;
   
for j=1:state
    for g=1:agent
ewsb(:,j,g)=wsb(:,:,g)*P(j,:)';
    end
end

if min(wsb)<-5000
c(:,:,:)=-5;   
countgrid(:,:)=0;
end   
if sum(sum(countgrid(i,j)))>0
 error('countgrid not zero')
end
PROB(:,:,1)=P;
STATE(:,:,1)=S;
WSB(:,:,:,Qbeta,Qrho,t)=wsb;
PHI(:,:,:,Qbeta,Qrho,t)=phi;
CONS(:,:,:,Qbeta,Qrho,t)=c;
GAP(Qbeta,Qrho,t)=gap;
XXKEEP(:,:,:,Qbeta,Qrho,t)=xkeep;
COUNTGRID(:,:,Qbeta,Qrho,t)=countgrid;
end % of iteration loop

PROBDEV3(1:state,1:state,1)=P;
STATEDEV3(1:state,1:agent,1)=S;
WDEV3(1:size(wsb,1),1:size(wsb,2),1:3,Qbeta,Qrho)=wsb;
EWDEV3(1:size(ewsb,1),1:size(ewsb,2),1:3,Qbeta,Qrho)=ewsb;
PHIDEV3(1:size(phisol,1),1:size(phisol,2),1:3,Qbeta,Qrho)=phisol;
CDEV3(1:size(csol,1),1:size(csol,2),1:3,Qbeta,Qrho)=csol;


save ('RESULTS3ALLNOAUT.mat', 'PROBDEV3','STATEDEV3','WDEV3','EWDEV3','PHIDEV3','CDEV3');


